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Summary: We report on a series of tests of Newtonian Lagrangian perturbation 
schemes using N-body simulations for various power-spectra with scale-independent 
indices in the range —3 to +1. The models have been evolved deeply into the non-linear 
regime of structure formation in order to probe the dynamical and statistical perfor- 
mance of the Lagrangian perturbation schemes (whose first-order solution contains as 
a subset the celebrated "Zel'dovich-approximation" (hereafter ZA). These tests reveal 
properties of the approximations at stages beyond the obvious validity of perturbation 
theory. Recently, another series of tests of different analytical and semi-numerical ap- 
proximations for large-scale structure was conducted with the result that ZA displays 
the best dynamical performance in comparison with the N-body simulations, if the ini- 
tial data were smoothed before evolving the model, i.e., a truncated form of ZA (TZA). 
We show in this Letter that the excellent performance of TZA can be further improved 
by going to second order in the Lagrangian perturbation approach. The truncated 
second-order Lagrangian scheme provides a useful improvement over TZA especially 
for negative power indices, which suggests it will be very useful for modelling standard 
scenarios such as "Cold-" , "Hot-" and "Mixed-Dark-Matter" . 
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1. Lagrangian perturbation theory put into perspective 



Zel'dovich (1970, 1973) proposed an approximation (hereafter ZA) by extrapolating the 
Eulerian linear theory of gravitational instability into the non-linear regime using the 
Lagrangian picture of continuum mechanics. He discussed interesting consequences of this 
approximation which is capable of describing shell-crossing singularities which, in this 
model, develop into highly anisotropic oblate "pancake" structures. Those structures can 
be fundamentally understood and classified in the framework of the Lagrange-singularity- 
theory, (Arnol'd et al. 1982, Shandarin and Zel'dovich 1989, Roshansky et al. 1994 and 
references therein). Zel'dovich's work initiated many applications making it one of the most 
cited articles in astronomy; it emerged as a standard tool to model principal elements 
of the large-scale structure, is used to initialize most N-body codes employed by the 
cosmology community, and forms the basis of more sophisticated approximations like the 
adhesion approximation (Gurbatov et al. 1989, Kofman et al. 1992). Zel'dovich's model 
can be derived by formulating the Euler-Poisson system in terms of Lagrangian coordinates 
and solving the Lagrangian evolution equations for the field of trajectories perturbatively 
(Buchert and Gotz 1987, Buchert 1989, 1992). It appears as a subclass of the irrotational 
Lagrangian first-order solution which covers substantial non-linearities in contrast to the 
Eulerian first-order solution. This explains the success of this approximation if applied to 
non-linear gravitational structure formation (compare Coles et al. 1993). 

The particular justification that ZA could be relevant to hierarchical clustering has 
developed slowly (Melott et al. 1983; Melott and Shandarin 1990; Kofman 1991; Little et 
al. 1991; Kofman et al. 1992; Coles et al. 1993; however see Peebles 1993). A general 
concensus is forming based around a unification of the former Soviet ("pancake") and 
Western (hierarchical clustering) theories. 

The Lagrangian theory of gravitational instability is now used in large-scale structure 
modelling much as the Eulerian theory of gravitational instability used to be. Numerous 
efforts concern the investigation and application of Lagrangian perturbation solutions up to 
the third order (Buchert 1989, 1992, Moutarde et al. 1991, Bouchet et al. 1992, Buchert 

1993, Buchert and Ehlers 1993, Gramann 1993, Giavalisco et al. 1993, Lachieze-Rey 
1993a,b, Buchert 1994, Juszkiewicz et al. 1994, Bernardeau 1994, Munshi and Starobinsky 

1994, Munshi et al. 1994), and, most recently, the investigation of general relativistic 
analogues (which are intrinsically Lagrangian in the eigensystem of the flow) (Matarrese 
et al. 1993, Kasai 1993, Croudace et al. 1994, Bertschinger and Jain 1994, Matarrese et 
al. 1994a,b, Berschinger and Hamilton 1994, Kofman and Pogosyan 1994, Salopek et al. 
1994). 
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2. Optimization of Lagrangian perturbation schemes 



Until recently, ZA was evolved only until shell-crossing, i.e., when singularities in the 
density field develop, at the epoch when the Eulerian representation of the basic dynamical 
equations breaks down. In principle, the Lagrangian representation of the flow allows 
following the evolution across caustics, where the flow field itself remains finite. This 
implies neglecting self-gravitating interaction of multi-stream systems developing inside 
caustics, however, secondary generations of shell-crosssings can be modelled as observed in 
N-body simulations by going to higher orders in the perturbation approach (Buchert and 
Ehlers 1993). Melott et al. 1994a (hereafter MPS) investigated the performance of a new 
approximation which requires truncation of high frequencies in the initial power-spectrum 
before evolving ZA (hereafter: TZA) taking the evolution of large-scale structure deeply 
into the non-linear regime. MPS found that filtering the initial data with a Gaussian at a 
scale close to but smaller than the non-linearity scale yields the best agreement with the 
density fields of the same (untruncated) initial data as evolved by an N-body code. 
The non-linearity scale k n i is defined by: 

fk nl {t) 

a 2 (t) / d 3 kV(k) = 1 , (1) 

Jo 

where k n i(t) is decreasing with time as successively larger scales enter the non-linear 
regime; a(t) is the scale factor of the homogeneous background (a(U) = 1), and V(k) 
denotes the initial power-spectrum taken to be a powerlaw with indices in the range —3 
to +1. 

"Best agreement" was defined in terms of an optimal scale k op t in /c-space at which 
the usual cross-correlation coefficient S between the resulting density fields attains its 
maximum: 

S := < > , (2) 

0-1(72 

where 5i,£ = 1,2 represent the density contrasts in the analytical and the numerical 
approximations, respectively, oz = a/< > — < Si > 2 is the standard deviation in a 
Gaussian random field; averages < ... > are taken over the entire distribution. We believe 
this is the most important statistical test, because it measures whether the approximation 
is moving mass to the right place, with an emphasis on dense regions. We also allow for 
small errors by calculating S for the two density arrays smoothed at a variety of smoothing 
lengths. 

For the Lagrangian perturbation schemes up to the third order which were used in our 
tests see (Buchert 1994). We conducted several tests: In the first step we studied "pancake 
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models", i.e., models which a priori have a truncated power-spectrum, in order to study 
principal effects of a second- and higher-order correction to ZA (for details see Buchert 
et al. 1994a). In the second step we analyzed the whole family of models with powerlaw- 
spectra —3, . . . , +1 by evolving them deeply into the non-linear regime (for details see 
Melott et al. 1994b). These "hierarchical models" have been evolved for expansion factors 
of 240 to 5100, depending on spectral index and k n i. 

Besides the cross-correlation coefficient (3) as a function of scale, we analyzed sev- 
eral statistics including the comparison of the evolved power-spectrum, the evolved r.m.s. 
values of the density contrast as a function of scale, the phase-angle accuracy achieved by 
the analytical models, and the evolved density distribution functions. 

For all these statistics and for all spectra studied with different filter types and filter 
scales, we always found improvement for the second-order scheme upon first-order (TZA), 
if the initial data are truncated with a Gaussian filter at a slightly larger scale k opt than 
the scale needed for the optimal TZA. 

In Fig. 1 we display slices of the final density fields for the spectrum with index 
n = — 1 as calculated by the N-body code and the first- and second-order Lagrangian 
approximations at the stage where the non-linearity scale has evolved to k n i =8kf, kf is 
the fundamental mode of the simulation box. 

In Fig. 2 we present the results of the most important statistical test which probes 
the dynamics of the models, i.e., the cross-correlation coefficient (3) for the whole family 
of hierarchical models at the same evolution stage. 

Considerably more details can be found in Melott et al. (1994b), where we show more 
statistical tests applied for the range of indices — 3 < n < 1. 
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3. Conclusions 



We summarize our main conclusions and list the advantages of going to second-order in 
perturbation theory for the purpose of modelling highly non-linear stages: 

• 1. The statistics which probe the gravitational dynamics of the models show im- 
provement due to second-order corrections. This success is found for a considerably higher 
non-linearity than expected from a perturbation approach. 

• 2. The improvement (although minor for much small-scale power) is robustbj going 
to later stages and to smaller scales. This holds for any spectrum and for any statistics 
analyzed. 

• 3. The CPU times on a CRAY YMP are for the first-order scheme 25 seconds, and 
for the second-order scheme 60 seconds; the corresponding CPU times on a CONVEX 
C220 are 2 and 5 minutes. Thus, even the second-order scheme is competitive with one 
step in a corresponding PM-type N-body simulation. 

• 4. The high speed as well as the fact that the second-order scheme is as easy to 
implement as the first-order scheme (directly from the initial data), render this model 
suitable for all areas of application where thus far ZA was used, e.g., the initialization of 
N-body codes. 

• 5. The second-order scheme predicts much faster collapse of first objects (treating 
also tidal effects) at times comparable to the collapse time in the widely used spherical 
"tophat" model (Moutarde et al 1991, Buchert and Ehlers 1993, Munshi et al. 1994). 
Thus, it is preferred for the treatment of ensembles of collapsing objects and for normal- 
ization purposes. 

• 6. Since the second-order corrections to TZA provide noticeable improvement of 
dynamical accuracy for initial data with negative sloped power-spectra, we expect that 
the truncated second-order scheme will be especially useful for the modelling of standard 
cosmogonies (like Hot-, Cold-, and Mixed-Dark-Matter). 

• 7. This modelling will be effective for large sample calculations, since in numerical 
realizations of 'fair' samples in excess of 300 h _1 Mpc, performed with the same resolution 
as the simulations reported here (128 3 particles on 128 3 meshs), the truncation scale is 
close to the Nyquist frequency of the N-body computing. Thus, shortcomings of the 
analytical schemes become negligible which puts them in an ideal position for the purpose 
of simulating the environment of galaxy formation down to scales where other physical 
effects start to affect models based on the description of self-gravity alone. Our method 
can be effective down to galaxy group mass scales (1O 13 M ), or better if we include biasing 
or go to epochs earlier than the present. Thus many things which have been studied by 
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N-body simulation can now be generated by approximation. The code for second-order is 
available on request from tob @ mpa-garching. mpg.de . 

• 8. The third-order scheme does not show the 'robustness' observed for second- 
and first-order. However, to draw definite conclusions the analytical solution of the third- 
order effect must be studied with reduced numerical uncertainties in its realization by 
Fast-Fourier-Transform, as pursued by Buchert et al. (1994b). 
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Figure Captions 

Figure 1: Thin slices (thickness L/128) of the density fields are displayed for the numerical 
(a), the optimally truncated first-order (b) and second-order (c) approximations for the evolution 
stage corresponding to k n i = 8k f, and for the single power-spectrum with index n = — 1. The 
grey-scale is logarithmic in order to emphasize the high-density regions. 

Figure 2: The cross-correlation coefficient S as a function of the standard deviation a p of 
the smoothed numerical simulation for the different power-spectra n = —3, —2, —1, 0, +1 (Figs. 
2a,b,c,d,e). The cross-correlation of the N-body with the optimally truncated first-order model 
is shown as a dotted line; with the optimally truncated second-order model a dashed line. 
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